#include "approximatedosfunctor.h"

const string ApproximateDoSFunctor::wavenumbersymbol = "WaveVector";
const string ApproximateDoSFunctor::statessymbol = "LanczosStates";
const string ApproximateDoSFunctor::etasymbol = "Eta";

ApproximateDoSFunctor::ApproximateDoSFunctor(System* s, const rule::list &opts) : hamFunc(s,opts), k(arma::fill::zeros), statetype(defaulttype), states(0), eta(0.001)
{
	const size_t rulesize = opts.size();
	for(size_t r = 0; r < rulesize; r++)
	{
		if(!opts[r].first.isSymbolType())
			continue;
		
		symbol sym = opts[r].first.toSymbol();
		variant var = opts[r].second;
		
		if(sym.name() == wavenumbersymbol)
			k = (var.isNumericListType()) ? var.toArmaDoubleVector3() : arma::vec3(arma::fill::zeros);
		else if(sym.name() == statessymbol && var.isNumericType())
		{
			statetype = userdefined;
			states = var.toLong();
		}
		else if(sym.name() == etasymbol)
			eta = (var.isNumericType()) ? var.toDouble() : 0.001;
	}
}

arma::cx_vec ApproximateDoSFunctor::generateRandomPhaseState(size_t n)
{
	arma::cx_vec result(n,arma::fill::zeros);
	result.set_imag(arma::randu<vec>(n));
	result = arma::exp(2*arma::datum::pi*result);
	return result/sqrt((double)n);
}

void ApproximateDoSFunctor::generateTriDiagHamElem(const arma::sp_cx_mat &ham, arma::cx_vec &alphas, arma::vec &betas, const arma::cx_vec &prevstate, arma::cx_vec &currentstate, long currentindex)
{
	arma::cx_vec H_currentstate = ham*currentstate;
	alphas(currentindex) = (currentstate.ht()*H_currentstate).eval()(0);
	if(currentindex == 0)
		currentstate = H_currentstate - alphas(currentindex)*currentstate;
	else
		currentstate = H_currentstate - alphas(currentindex)*currentstate - betas(currentindex-1)*prevstate;
	betas(currentindex) = arma::norm(currentstate,2);
	currentstate = currentstate/betas(currentindex);
}

complex ApproximateDoSFunctor::generateGreensResult(const arma::cx_vec &alphas, const arma::vec &betas, const complex &energy, long currentindex, complex &lastGreen)
{
	if(currentindex == betas.size() -1)
	{
		lastGreen = (energy - alphas(currentindex) - sqrt(pow(energy-alphas(currentindex),2) - 4. * pow(betas(currentindex-1),2) ) )/( 2. * pow(betas(currentindex-1),2) );
		if(lastGreen.imag() > 0)
			return lastGreen;
	
		return (lastGreen = (energy - alphas(currentindex) + sqrt(pow(energy-alphas(currentindex),2) - 4. * pow(betas(currentindex-1),2) ) )/( 2. * pow(betas(currentindex-1),2) ));
	}
	else
		return (lastGreen = 1./(energy - alphas(currentindex) - pow(betas(currentindex),2)*lastGreen));
}

complex ApproximateDoSFunctor::getGreensFunctionValue(int id, const complex &energy)
{
	arma::sp_cx_mat ham = hamFunc.getTransferIntegral(id,k);

	size_t numberofstates = (statetype == defaulttype) ? 1000 : states;

	alphas.set_size(numberofstates);
	betas.set_size(numberofstates);
	arma::cx_vec firststate = generateRandomPhaseState(ham.n_rows);
	arma::cx_vec prevstate;
	arma::cx_vec tmp;

	for(long i = 0; i < numberofstates; i++)
	{
		prevstate = tmp;
		tmp = firststate;
		generateTriDiagHamElem(ham,alphas,betas,prevstate,firststate,i);
	}
	
	complex currentval = 0.;
	for(long i = numberofstates-1; i >= 0; i--)
		generateGreensResult(alphas, betas, energy, i, currentval);

	return currentval;
}

double ApproximateDoSFunctor::getDoSValue(int id, double energy)
{
	complex Z = complex(energy,eta);
	return -getGreensFunctionValue(id, Z).imag()/arma::datum::pi;
}

arma::mat ApproximateDoSFunctor::getDoS(int id, numeric::range energies)
{
	arma::sp_cx_mat ham = hamFunc.getTransferIntegral(id,k);

	size_t numberofstates = (statetype == defaulttype) ? 1000 : states;

	alphas.set_size(numberofstates);
	betas.set_size(numberofstates);
	const arma::cx_vec firststate = generateRandomPhaseState(ham.n_rows);

	arma::cx_vec prevstate,tmp,currentstate;
	currentstate = firststate;

	for(long i = 0; i < numberofstates; i++)
	{
		prevstate = tmp;
		tmp = currentstate;
		generateTriDiagHamElem(ham,alphas,betas,prevstate,currentstate,i);
	}
	
	const size_t size = energies.size();
	arma::mat result(2,size);

	_Cilk_for(size_t e = 0; e < size; e++)
	{
		complex energy = complex((double)energies.begin + e*(double)energies.step,eta);
		
		complex currentval = 0.;
		for(long i = numberofstates-1; i >= 0; i--)
			generateGreensResult(alphas, betas, energy, i, currentval);

		result(0,e) = energy.real();
		result(1,e) = -currentval.imag()/arma::datum::pi;
	}

	return result;
}